home *** CD-ROM | disk | FTP | other *** search
/ MacFormat España 15 / macformat_15.iso / Shareware Internet / Desarrolladores / gray image 2.1 / median_filter.cc < prev    next >
Text File  |  1995-05-31  |  9KB  |  297 lines
  1. // This may look like C code, but it is really -*- C++ -*-
  2. /*
  3.  ************************************************************************
  4.  *
  5.  *               Grayscale Image
  6.  *
  7.  *        Perform median filtration of an image
  8.  *
  9.  * Implementations of several flavors of median filtration: by rows
  10.  * only, filter only columns, filter rows and then filter columns. The
  11.  * latter is equivalent to a 2D median filtration with a diamond-shaped
  12.  * window. 2D filtration with a square window fits the existing framework,
  13.  * too, but not yet implemented.
  14.  *
  15.  * When the filtration window sticks out of the image, we extrapolate
  16.  * the corresponding pixels. Incidentally, it means that the border
  17.  * pixel(s) are never affected by the filtration. Indeed, say, for a
  18.  * 3-point row filtration
  19.  *    new_pixel[i] = median_of(pixel[i-1],pixel[i],pixel[i+1])
  20.  * if we assume that pixel[-1] = pixel[0], then
  21.  *    median_of(pixel[0],pixel[0],pixel[1]) is always pixel[0]
  22.  * (majority wins). The same holds for other types of filtration.
  23.  *
  24.  * Square window 2D filtration:
  25.  *    For every pixel of the source image x[i,j] compute the
  26.  *    pixel y[i,j] of the filtrated image according to the
  27.  *    rule
  28.  *    y[i,j] = median_of( x[i+p,j+l] ), p,l = -ws...0...ws
  29.  *
  30.  *    where 2*ws+1 is the window size. The median of the sequence
  31.  *    is the middle term of the ordered sequence (in our case the
  32.  *    sequence always contains the odd number of terms).
  33.  * This isn't implemented yet.
  34.  *
  35.  * $Id: median_filter.cc,v 2.5 1995/03/14 22:11:03 oleg Exp oleg $
  36.  *
  37.  ************************************************************************
  38.  */
  39.  
  40. #include "filter.h"
  41. #include "minmax.h"
  42.  
  43. /*
  44.  *------------------------------------------------------------------------
  45.  *        Simple 3-point median filtration
  46.  */
  47.  
  48.                 // Find a median of three numbers
  49. static inline GRAY_SIGNED median_of(const GRAY_SIGNED a1,
  50.                     const GRAY_SIGNED a2,
  51.                     const GRAY_SIGNED a3)
  52. {
  53.   if( a3 >= a1 )
  54.     if( a2 >= a3 )
  55.       return a3;            // order is a1, a3, a2
  56.     else
  57.       if( a2 > a1 )
  58.     return a2;            // order is a1, a2, a3
  59.       else
  60.     return a1;            // order is a2, a1, a3
  61.   else
  62.     if( a2 >= a1 )
  63.       return a1;            // order is a3, a1, a2
  64.     else
  65.       if( a2 > a3 )
  66.     return a2;            // order is a3, a2, a1
  67.       else
  68.     return a3;            // order is a2, a3, a1
  69. }
  70.  
  71.  
  72.                 // 3-point row-wise filtration
  73. IMAGE& FilterIt::median_3_row(void)
  74. {
  75.   if( image.ncols < 2 )
  76.     image.info(),
  77.     _error("we need at least 3 columns for 3-point median row filtration");
  78.  
  79.   register GRAY_SIGNED * pp = (GRAY_SIGNED *)image.pixels;
  80.   register GRAY_SIGNED prev_pixel;    // in the row
  81.  
  82.   for(GRAY_SIGNED * row_end = (GRAY_SIGNED *)image.pixels + image.ncols;
  83.       row_end <= (GRAY_SIGNED *)image.pixels + image.npixels;
  84.       row_end += image.ncols)
  85.   {
  86.     prev_pixel = *pp++;            // At the beginning, pixel[-1]=pixel[0]
  87.                     // and the first pixel never changes
  88.     while(pp < row_end-1)
  89.     {
  90.       GRAY_SIGNED new_pixel = median_of(prev_pixel,pp[0],pp[1]);
  91.       prev_pixel = pp[0];
  92.       *pp++ = new_pixel;
  93.     }
  94.     pp++;                // if we assume p[ncols] = p[ncols-1]
  95.                     // then the last pixel is always the
  96.                     // same
  97.   }
  98.   assert( pp == (GRAY_SIGNED *)image.pixels + image.npixels );
  99.   return image;
  100. }
  101.  
  102.                 // 3-point col-wise filtration
  103. IMAGE& FilterIt::median_3_col(void)
  104. {
  105.   if( image.nrows < 2 )
  106.     image.info(),
  107.     _error("we need at least 3 rows for 3-point median col filtration");
  108.  
  109.   register const int ncols = image.ncols;    // caching
  110.   const int fl_jump = image.npixels - ncols;    // between 1. and last rows
  111.   register GRAY_SIGNED * pp = (GRAY_SIGNED *)image.pixels;
  112.   register GRAY_SIGNED prev_pixel;    // in the col
  113.  
  114.   GRAY_SIGNED * last_row = (GRAY_SIGNED *)image.pixels + fl_jump;
  115.   for(; last_row < (GRAY_SIGNED *)image.pixels + image.npixels; last_row++) 
  116.   {
  117.     prev_pixel = *pp, pp += ncols;    // At the beginning, pixel[-1]=pixel[0]
  118.                     // and the first pixel always stays
  119.     for(; pp < last_row; pp += ncols)
  120.     {
  121.       GRAY_SIGNED new_pixel = median_of(prev_pixel,pp[0],pp[ncols]);
  122.       prev_pixel = pp[0];
  123.       pp[0] = new_pixel;
  124.     }
  125.     pp -= fl_jump-1;            // if we assume p[nrows] = p[nrows-1]
  126.                     // then the last pixel is always the
  127.                     // same
  128.   }
  129.   assert( pp == last_row - fl_jump );
  130.   return image;
  131. }
  132.  
  133. /*
  134.  *------------------------------------------------------------------------
  135.  *            Simple 5-point median filtration
  136.  *
  137.  * Note about the algorithm:
  138.  * we can always reduce finding a median of 5 numbers to finding a median
  139.  * of three numbers by noting that max and min of four numbers never
  140.  * count. Indeed, if mx = max(a1,a2,a4,a5) and a3 > mx, then
  141.  * the order is "? ? ? mx a3" and mx won't figure in in the median, so
  142.  * the median is the same as of "? ? ? a3". If a3 < mx, then the order
  143.  * is "? ? ? a3 mx" and mx won't count even more. The same holds for the
  144.  * min of four numbers.
  145.  */
  146.                 // Find a median of 5 numbers
  147. static inline GRAY_SIGNED median_of(const GRAY_SIGNED a1,
  148.                     const GRAY_SIGNED a2,
  149.                     const GRAY_SIGNED a3,
  150.                     const GRAY_SIGNED a4,
  151.                     const GRAY_SIGNED a5)
  152. {
  153.   register GRAY_SIGNED mn1 = a1, mx1 = a2;
  154.   if( a1 > a2 )
  155.     mn1 = a2, mx1 = a1;
  156.   register GRAY_SIGNED mn2 = a4, mx2 = a5;
  157.   if( a4 > a5 )
  158.     mn2 = a5, mx2 = a4;
  159.   return median_of(min(mx1,mx2),a3,max(mn1,mn2));
  160. }
  161.  
  162.                 // 5-point row-wise filtration
  163. IMAGE& FilterIt::median_5_row(void)
  164. {
  165.   if( image.ncols < 4 )
  166.     image.info(),
  167.     _error("we need at least 5 columns for 5-point median row filtration");
  168.  
  169.   register GRAY_SIGNED * pp = (GRAY_SIGNED *)image.pixels;
  170.   register GRAY_SIGNED prev_pixel1;    // pixel[*,j-1]
  171.   register GRAY_SIGNED prev_pixel2;    // pixel[*,j-2]
  172.  
  173.   for(GRAY_SIGNED * row_end = (GRAY_SIGNED *)image.pixels + image.ncols;
  174.       row_end <= (GRAY_SIGNED *)image.pixels + image.npixels;
  175.       row_end += image.ncols)
  176.   {
  177.     prev_pixel2 = prev_pixel1 = *pp++;    // pixel[-2] = pixel[-1] = pixel[0]
  178.                     // and the first pixel always stays
  179.     while(pp < row_end-2)
  180.     {
  181.       GRAY_SIGNED new_pixel = 
  182.     median_of(prev_pixel2,prev_pixel1,pp[0],pp[1],pp[2]);
  183.       prev_pixel2 = prev_pixel1;
  184.       prev_pixel1 = pp[0];
  185.       *pp++ = new_pixel;
  186.     }
  187.     *pp = median_of(prev_pixel2,prev_pixel1,pp[0],pp[1],pp[1]);
  188.     pp += 2;                // pp[ncols-1] is unchanged
  189.   }
  190.   assert( pp == (GRAY_SIGNED *)image.pixels + image.npixels );
  191.   return image;
  192. }
  193.  
  194.                 // 5-point col-wise filtration
  195. IMAGE& FilterIt::median_5_col(void)
  196. {
  197.   if( image.nrows < 4 )
  198.     image.info(),
  199.     _error("we need at least 5 rows for 5-point median col filtration");
  200.  
  201.   register const int ncols = image.ncols;    // caching
  202.   const int fl_jump = image.npixels - 2*ncols;    // between 1. and last b1 rows
  203.   register GRAY_SIGNED * pp = (GRAY_SIGNED *)image.pixels;
  204.   register GRAY_SIGNED prev_pixel1;    // pixel[i-1,*]
  205.   register GRAY_SIGNED prev_pixel2;    // pixel[i-2,*]
  206.   register GRAY_SIGNED next_pixel1;    // pixel[i+1,*]
  207.   register GRAY_SIGNED next_pixel2;    // pixel[i+2,*]
  208.  
  209.   GRAY_SIGNED * last_b1_row = (GRAY_SIGNED *)image.pixels + fl_jump;
  210.   for(; last_b1_row < (GRAY_SIGNED *)image.pixels + image.npixels - ncols;
  211.       last_b1_row++) 
  212.   {
  213.     prev_pixel1 = prev_pixel2 = *pp, pp += ncols; // skip over pixel[0,*]
  214.                           // which is unaffected
  215.     next_pixel1 = pp[ncols];
  216.     for(; pp < last_b1_row; pp += ncols)
  217.     {
  218.       next_pixel2 = pp[2*ncols];
  219.       GRAY_SIGNED new_pixel = 
  220.     median_of(prev_pixel2,prev_pixel1,pp[0],next_pixel1,next_pixel2);
  221.       prev_pixel2 = prev_pixel1;
  222.       prev_pixel1 = pp[0];
  223.       next_pixel1 = next_pixel2;
  224.       pp[0] = new_pixel;
  225.     }
  226.     pp[0] = median_of(prev_pixel2,prev_pixel1,pp[0],next_pixel1,next_pixel1);
  227.     pp -= fl_jump-1;            // skip over pixel[nrows-1]
  228.   }
  229.   assert( pp == last_b1_row - fl_jump );
  230.   return image;
  231. }
  232.  
  233. /*
  234.  *------------------------------------------------------------------------
  235.  *            Median Filtration drivers
  236.  * pick up and execute the most efficient method given input parameters
  237.  */
  238.  
  239. IMAGE& FilterIt::median(const Direction how, const int win_size)
  240. {
  241.   image.is_valid();
  242.  
  243.   switch(win_size)
  244.   {
  245.     case 3:
  246.          switch(how)
  247.      {
  248.        case RowsAndColumns:
  249.             return median_3_row(), median_3_col();
  250.        case Columns:
  251.             return median_3_col();
  252.        case Rows:
  253.             return median_3_row();
  254.      };
  255.      break;
  256.  
  257.     case 5:
  258.          switch(how)
  259.      {
  260.        case RowsAndColumns:
  261.             return median_5_row(), median_5_col();
  262.        case Columns:
  263.             return median_5_col();
  264.        case Rows:
  265.             return median_5_row();
  266.      };
  267.      break;
  268.  
  269.     default:
  270.          _error("Sorry, not implemented yet. "
  271.         "FYI: implemented windows: 3 and 5, and they are fast");
  272.   }
  273.   _error("There is no way we can reach here...");
  274.   return image;                // just to make the compiler happy
  275. }
  276.  
  277. #if 0
  278. /*
  279.  *------------------------------------------------------------------------
  280.  *        Finding a median of the sequence of GRAYs
  281.  * 
  282.  * Algorithm
  283.  *    Described in the book "Numerical Recepies".
  284.  * The algorithm does not require sorting the sequence. The main idea is
  285.  * that if the median of the sequence is a "middle" term of the sequence,
  286.  * the following relation should hold
  287.  *
  288.  *    SUM{ (xmed - x[i]) / |xmed - x[i]| } = 0
  289.  *
  290.  */
  291.  
  292. static GRAY find_median(GRAY seq[], const int no_elems)
  293. {
  294. }
  295.  
  296. #endif
  297.